Problem Statement¶

You are interested in studying the factors that influence the likelihood of heart disease among patients.

You have a dataset of 303 patients, each with 14 variables: age, sex, chest pain type, resting blood pressure, serum cholesterol, fasting blood sugar, resting electrocardiographic results, maximum heart rate achieved, exercise induced angina, oldpeak, slope, number of major vessels, thalassemia, and diagnosis of heart disease.

You want to use Bayesian logistic regression to model the probability of heart disease (the outcome variable) as a function of some or all of the other variables (the predictor variables).

You also want to compare different models and assess their fit and predictive performance.

Bayesian Workflow¶

For this project, I will try to follow this workflow:

1. Data exploration: Explore the data using descriptive statistics and visualizations to get a sense of the distribution, range, and correlation of the variables. Identify any outliers, missing values, or potential errors in the data. Transform or standardize the variables if needed.

2. Model specification: Specify a probabilistic model that relates the outcome variable to the predictor variables using a logistic regression equation. Choose appropriate priors for the model parameters, such as normal, student-t, or Cauchy distributions. You can use the brms package in Julia to define and fit Bayesian models using a formula syntax similar to lme4. However, try to use Turing.jl

3. Model fitting: Fit the model using a sampling algorithm such as Hamiltonian Monte Carlo (HMC) or No-U-Turn Sampler (NUTS). You can use the DynamicHMC or Turing.jl package in Julia to implement these algorithms. Check the convergence and mixing of the chains using diagnostics such as trace plots, autocorrelation plots, effective sample size, and potential scale reduction factor. You can use the MCMCDiagnostics or the included diagnostics features in Turing.jl package in Julia to compute these diagnostics.

4. Model checking: Check the fit and validity of the model using posterior predictive checks, residual analysis, and sensitivity analysis. You can use the PPCheck package in Julia to perform posterior predictive checks, which compare the observed data to data simulated from the posterior predictive distribution. You can use the BayesianRidgeRegression package in Julia to perform residual analysis, which plots the residuals against the fitted values and the predictor variables. You can use the Sensitivity package in Julia to perform sensitivity analysis, which measures how the posterior distribution changes with respect to the prior distribution or the likelihood function.

Data Exploration¶

After "collecting" the data, we may import it and arrange it so we can use it further.

The data set can be found in this Kaggle link.

In the above data frame, the attributes are as follows:

Complete attribute documentation:

1. age: age in years
2. sex: sex (1 = male; 0 = female)
3. cp: chest pain type
	- Value 1: typical angina
	- Value 2: atypical angina
	- Value 3: non-anginal pain
	- Value 4: asymptomatic
4. trestbps: resting blood pressure (in mm Hg on admission to the
hospital)
5. chol: serum cholestoral in mg/dl
6.fbs: fasting blood sugar > 120 mg/dl (1 = true; 0 = false)
7. restecg: resting electrocardiographic results
	- Value 0: normal
	- Value 1: having ST-T wave abnormality (T wave inversions and/or ST elevation or depression of > 0.05 mV)
	- Value 2: showing probable or definite left ventricular hypertrophy by Estes' criteria
8. thalach: maximum heart rate achieved
9. exang: exercise induced angina (1 = yes; 0 = no)
10. oldpeak: ST depression induced by exercise relative to rest
11. slope: the slope of the peak exercise ST segment
	- Value 1: upsloping
	- Value 2: flat
	- Value 3: downsloping
12. ca: number of major vessels (0-3) colored by flourosopy (for calcification of vessels)
13. thal: results of nuclear stress test (3 = normal; 6 = fixed defect; 7 = reversable defect)
14. num: target variable representing diagnosis of heart disease (angiographic disease status) in any major vessel
	- Value 0: < 50% diameter narrowing
	- Value 1: > 50% diameter narrowing

Data Interpretation¶

After collecting the data, it has been imported as a Data Frame. Now, to understand what we will do with this exercise, we need to analyze the data by means of Bayesian Logistic Regression.

With this type of analysis, we can make predictions on (typically) binary outcomes, based on a set of parameters. In this particular case, we are interested in predicting whether a patient will have heart disease based on a set of parameters such as age, chest pain, blood pressure, etc.

In terms of the data available, we have a set of 303 observations (303 patients) whose symptoms and circumstances have been recorded, and the outcome is the heart disease diagnosis. To simplify things, this data set has a binary outcome, i.e. heart disease present/not present.

Additionally, this study is divided in two parts: first, I will set up the logistic regression model to include only one predictor, i.e., age. Afterwards, an analysis will be performed including two or more predictors.

Model Specification¶

In this stage of the workflow, we will specify the Bayesian model and then use Turing.jl to program it in Julia.

The model I will use for this analysis is a Bayesian Logistic Regression model, which relates the probability of heart disease to a linear combination of the predictor variables, using a logistic function. The model can be written as:

where $y_i$ is the outcome for the i-th patient, $p_i$ is the probability of heart disease for the i-th patient, $\eta_i$ is the linear predictor for the i-th patient, $\alpha$ and $\beta_j$ are the intercept and coefficient for the j-th predictor variable, respectively, and $x_{ij}$ is the value of the j-th predictor variable for the i-th patient.

The assumptions that I am making are:

1. The outcome variable follows a Bernoulli distribution, i.e. $y_i \sim Bernoulli(p_i)$, which is appropriate for binary outcomes
2. The predictor variables are linearly related to the log-odds of the outcome variable, i.e. $\log(\frac{p}{1-p})$ which is a common assumption for logistic regression models
3. The prior distributions for the model parameters are uniform, which are weakly informative and reflect my prior beliefs about the plausible range of the parameters

Regarding point (2):

That statement means that the log-odds of the outcome variable (the log of the odds ratio) can be expressed as a linear function of the predictor variables. Mathematically, this can be written as:

where $p$ is the probability of the outcome variable being 1, $x_1, x_2, \ldots, x_k$ are the predictor variables, and $\alpha, \beta_1, \beta_2, \ldots, \beta_k$ are the coefficients (parameters).

This assumption implies that the effect of each predictor variable on the log-odds of the outcome variable is contant, regardless of the values of the other predictor variables. It also implies that the relationship between the predictor variables and the probability of the outcome variable is non-linear, as the probability is obtained by applying the inverse of the log-odds function, which is the logistic function:

The logistic function is an S-shaped curve that maps any real number to a value between 0 and 1. It has the property that as the linear predictor increases, the probability approaches 1, and as the linear predictor decreases, the probability approaches 0.

Model Specification Using Turing.jl¶

Crank up the Bayes!¶

Run the model using sample(model, sampler, opt_argument, samples, chains)

NOTE: The above routine employs the MCMCThreads() method to sample multiple chains. However, to implement this, one needs to change the environment variables for the number of threads Julia can use. These two discussions might shed some light as to how to achieve this:

1. https://docs.julialang.org/en/v1/manual/multi-threading/#man-multithreading
2. https://discourse.julialang.org/t/julia-num-threads-in-vs-code-windows-10-wsl/28794

Of course, if you don't want to bother, then just change the last two functional lines in the cell above so that they read:

# set number of chains - comment this out:
# num_chains = 8

# crank up the Bayes! - delete MCMCThreads() and num_chains chain = sample(model, sampler, samples)

Plot the MCMC Diagnostics¶

Get the Summary Statistics¶

Plot and Interpret the Results¶

Ok, how do we interpret the results from a Bayesian approach? Let's start by plotting the results. This will help us understand not only the results, but really grasp the power of a Bayesian model in action.

From a frequentist or a machine learning approach, we would expect to find a function that models the data the best possible way, i.e. fit a model. If we were to visualize it, we would see one single sigmoid curve trying its best to explain the data.

How about this chart here, though? This chart is a collection of possible outcomes given that the parameters $\alpha$ and $\beta$ in this case, are modeled as random variables with some probability distribution. Therefore, there is an uncertainty associated with them. This uncertainty is naturally propagated onto the sigmoid function. Therefore, there is also an uncertainty associated with that sigmoid curve that we are trying to model.

Again, below we can see a collection of possible outcomes given the parameter sample space. There is a darker region where most sigmoid functions turned out, and these tend to be the most probable sigmoid functions, or, in other words, these sigmoid functions are the most probable functions that could fit the data, considering the distributions of the parameters too!

Making Predictions¶

So why go through all this trouble, you might be asking. Well, one of the reasons we want to use probabilistic models is, first, to make predictions. But I would go further than that: these models are useful when making informed decisions. Let's try this out.

Let's make predictions for different arbitrary ages (50, 60, 70, 80, 20):

Interpreting the predictions¶

The last operations make predictions of heart diseased based on age only. What the predictions mean is that, given the data, the probability distribution of an individual of age 50 to have heart disease has a mean of 0.379, and a standard deviation of 0.485 (this is highly uncertain, by the way).

Similarly, a 20-year-old individual has a probability with a mean of 0.13 and standard deviation of 0.336 of having heart disease.

These statistics are extremely powerful when you are trying to make decisions, such as when diagnosing Heart Disease. It stands to reason that, if you were a physician, you want to know what your model says might be wrong (or not) with your patient, but you also want to know how much you can trust that prediction.

If your model classifies Patient X as having heart disease, you would probably want to know how sure you are of this. And this certainty comes partially from... you guessed it: your priors and the data.

In the plot below, we can see the where the predictions lie. Note that these probabilities are on a continuum given by the sigmoid function. But we want our final decision to be a yes or a no. To do that, we need to set a decision threshold.

We will do that at the end of the next section.

Model Specification Using Multiple Predictors¶

Some Data Cleaning¶

In this part, I am using the Turing.jl documentation tutorial found in https://turinglang.org/dev/tutorials/02-logistic-regression/.

In the tutorial, they quite rightly incorporate a train/test split, and data normalization, which is the recommended practice. I didn't do it in the first part of this tutorial to keep things simple!

Here is how they handle the split and the data normalization using MLUtils.

function split_data(df, target; at=0.70)
    shuffled = shuffleobs(df)
    return trainset, testset = stratifiedobs(row -> row[target], shuffled; p=at)
end

features = [:StudentNum, :Balance, :Income] numerics = [:Balance, :Income] target = :DefaultNum
trainset, testset = split_data(data, target; at=0.05) for feature in numerics μ, σ = rescale!(trainset[!, feature]; obsdim=1) rescale!(testset[!, feature], μ, σ; obsdim=1) end
# Turing requires data in matrix form, not dataframe train = Matrix(trainset[:, features]) test = Matrix(testset[:, features]) train_label = trainset[:, target] test_label = testset[:, target];

Inference¶

Now that our data is formatted, we can perform our Bayesian logistic regression with multiple predictors: using chest pain (cp), age (age), resting bloodpressure (tresttbps) and cholesterol (chol) levels.

Now we build the model and create the chain:

Plot the MCMC Diagnostics¶

Summary Statistics¶

Thank you!¶

And that concludes this little tutorial showcasing the power of a Bayesian model and the fun of using Julia. Thank you for stopping by!

Victor Flores